Step by step Time Homogeneous Markov Model in heemod: PCI vs CABG

1 Model description

We compare Percutaneous Coronary Intervention (PCI) vs Coronary Artery Bypass Grafting (CABG) over 10 annual cycles in a time-homogeneous cohort Markov model.

2 Transition probabilities - how patients move between states

Idea: Each row of the transition matrix is “from state”; each column is “to state”. Rows must sum to 1. Death’s row is [0,0,0,1].

Function you’ll see: define_transition() builds a transition object that strategies will use.

library(heemod)

# State names
nm <- c("Stable","MI","Reintervention","Death")

# --- PCI transition matrix (rows = from, columns = to) ---
mat_pci <- define_transition(
  
    0.70, 0.15, 0.10, 0.05,  # from Stable
    0.50, 0.00, 0.20, 0.30,  # from MI
    0.80, 0.00, 0.00, 0.20,  # from Reintervention
    0.00, 0.00, 0.00, 1.00   # from Death (absorbing)
  ,
  state_names = nm
)


# --- CABG transition matrix ---
mat_cabg <- define_transition(
    0.80, 0.08, 0.07, 0.05,  # from Stable
    0.60, 0.00, 0.10, 0.30,  # from MI
    0.80, 0.00, 0.00, 0.20,  # from Reintervention
    0.00, 0.00, 0.00, 1.00   # from Death
  ,
  state_names = nm
)

Common pitfalls: Rows not summing to 1; mismatched state names between matrices and strategies; Death row not absorbing.

3 Plot model

3.1 PCI

plot(mat_pci)

3.2 CABG

plot(mat_cabg)

4 State values

Idea: Each state accrues cost and QALY for the time spent there in a cycle.

Function you’ll see: define_state() sets these values. Here we discount inside states using discount(rate). All states must expose the same value names.

# Economic parameters
cost_discount_rate <- 0.03
qaly_discount_rate <- 0.03

cost_stable <-  50000
cost_mi     <- 20000
cost_reint  <- 150000
cost_death  <- 0

util_stable <- 0.85
util_mi     <- 0.60
util_reint  <- 0.70
util_death  <- 0.00

# Define states (cost_total = state cost + strategy-specific upfront (only in cycle 1))
# Define states: discount inside states (no cycle logic)
state_Stable <- define_state(
  cost_total = discount(cost_stable, cost_discount_rate),
  qaly       = discount(util_stable, qaly_discount_rate)
)

state_MI <- define_state(
  cost_total = discount(cost_mi, cost_discount_rate),
  qaly       = discount(util_mi, qaly_discount_rate)
)

state_Reintervention <- define_state(
  cost_total = discount(cost_reint, cost_discount_rate),
  qaly       = discount(util_reint, qaly_discount_rate)
)

state_Death <- define_state(
  cost_total = 0,
  qaly       = discount(util_death, qaly_discount_rate)
)

Common pitfalls: Inconsistent value names across states; double discounting (don’t also pass discount= to run_model() if you discount here).

One-time upfront costs — cleanly at baseline

Idea: Charge the procedure once at the start (cycle 0), not every cycle.

Function you’ll see: define_starting_values() inside define_strategy() adds one-time values at baseline. The names must match what we pass to run_model(cost=..., effect=...).

# One-time costs at baseline
upfront_pci  <- 200000
upfront_cabg <- 300000

5 Strategy definitions — assemble transitions, states, and upfronts

Function you’ll see: define_strategy(transition=..., <state>=..., starting_values=...)

strat_pci <- define_strategy(
  transition = mat_pci,
  Stable = state_Stable,
  MI = state_MI,
  Reintervention = state_Reintervention,
  Death = state_Death,
  starting_values = define_starting_values(
    cost_total = upfront_pci,  # one-time baseline cost
    qaly       = 0             # no baseline QALY
  )
)

strat_cabg <- define_strategy(
  transition = mat_cabg,
  Stable = state_Stable,
  MI = state_MI,
  Reintervention = state_Reintervention,
  Death = state_Death,
  starting_values = define_starting_values(
    cost_total = upfront_cabg,  # one-time baseline cost
    qaly       = 0             # no baseline QALY
  )
)

Common pitfalls: Strategy states must match transition state names and order; all states must expose the same value names.

6 Running the model— simulate the cohort

Function you’ll see: run_model() executes the Markov model across cycles.

  • We choose method = "end" (values counted at cycle end).
  • Because we discounted inside states, we do not pass discount= here.
res <- run_model(
  pci  = strat_pci,
  cabg = strat_cabg,
  cycles = 10,
  init = c(Stable = 1000, MI = 0, Reintervention = 0, Death = 0),
  method = "beginning",
  cost = cost_total,
  effect = qaly
)

Alternative pattern: Put raw values in states and pass discount = c(cost=0.05, effect=0.05) here (but then don’t use discount() in states).

7 Read results — totals, ICER, and NMB

Functions you’ll see: summary() (version‑dependent). The totals already include the one‑time starting values we set earlier.

wtp <- c(150000, 300000)
summary(res, threshold = wtp)
2 strategies run for 10 cycles.

Initial state counts:

Stable = 1000
MI = 0
Reintervention = 0
Death = 0

Counting method: 'beginning'.

Values:

     cost_total     qaly
pci   543798010 4957.920
cabg  653200098 5367.794

Net monetary benefit difference:

       150000    3e+05
pci  47921.02     0.00
cabg     0.00 13560.05

Efficiency frontier:

pci -> cabg

Differences:

     Cost Diff. Effect Diff.     ICER Ref.
cabg   109402.1    0.4098738 266916.5  pci

Formulas to remember:
- ICER = (Cost_CABG − Cost_PCI) / (QALY_CABG − QALY_PCI)
- NMB = QALY * WTP − Cost

8 Visualise — see how the cohort and values evolve

Function you’ll see: plot(res, ...) with built‑in panels.

8.1 Patient Counts by Strategy

library(ggplot2)
plot(res, type = "counts")+theme_minimal()

8.2 Patient Counts by State

plot(res, type = "counts", panel ="by_state", free_y = TRUE)+
  theme_minimal()

8.3 Costs by Strategy

plot(res, type = "values", panel = "by_value",
     free_y = TRUE) +
  scale_color_brewer(
    name = "Strategy",
    palette = "Set1"
  )+
  theme_minimal()

9 Trace tables

9.1 Population over cycles

population_over_cycles <- get_counts(res)

t_population <- population_over_cycles |> 
  tidyr::pivot_wider(names_from = .strategy_names, values_from = count) |>
  dplyr::arrange(model_time)

library(DT)
DT::datatable(t_population,
              extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
                pageLength = 5
              ))

9.2 Costs and QALYs over cycles

costs_qalys_over_cycles <- get_values(res)
t_costs_qalys <- costs_qalys_over_cycles |> 
  tidyr::pivot_wider(names_from = .strategy_names, values_from = value) |>
  dplyr::arrange(model_time, value_names)
DT::datatable(t_costs_qalys,
              extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
                pageLength = 5
              ))